% 在固定 mu=3.5, z0=-1 条件下，通过改变耦合强度 k，
% 在 (x0,y0) 初始平面上扫描，绘制吸引域分布(basins of attraction)
%
% ​蓝色区域​（x≈0.2）：系统收敛到低值周期轨道
% ​红色区域​（x≈0.8）：系统收敛到高值周期轨道
% ​过渡色区域​（黄绿色）：混沌或复杂多周期状态
%
% 未对最终轨迹进行判定。


clc; clear; close all;

% 参数固定
mu = 3.5;       % 固定 mu
z0 = -1;        % 固定 z0
k_list = [0.15, 0.2, 0.25, 0.3];  % 不同耦合强度 (与论文 Fig.5 对应)
labels = {'(a) k=0.15', '(b) k=0.2', '(c) k=0.25', '(d) k=0.3'};

% 初始条件平面范围
Nx = 100; Ny = 100;           % 网格分辨率
x0_range = linspace(0, 1, Nx);
y0_range = linspace(0, 1, Ny);

% 迭代步数设置
nTransient = 200;  % 丢弃瞬态步数
nIter = 300;       % 正式迭代步数 (最终颜色根据最后一次迭代的 x 值)


for idxK = 1:4
    k = k_list(idxK);

    % 为每个 k 创建一个 Nx×Ny 的矩阵存储"最终 x"
    basinData = zeros(Ny, Nx);

    % 在 (x0,y0) 平面做网格扫描
    for ix = 1:Nx
        for iy = 1:Ny
            % (a) 初始化
            state = [x0_range(ix), y0_range(iy), z0];

            % (b) 丢弃瞬态
            for n = 1:nTransient
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
            end

            % (c) 正式迭代 nIter 步, 取最后一次的 x 值
            for n = 1:nIter
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
            end
            final_x = state(1);

            % (d) 收集 "最终 x 值" 
            basinData(iy, ix) = final_x;
        end
    end

    subplot(2,2,idxK);
    % 这里使用 imagesc 映射 "最终 x 值" => 颜色
    % x0_range 对应横轴, y0_range 对应纵轴
    imagesc(x0_range, y0_range, basinData);
    set(gca, 'YDir','normal');  % 使 y 轴向上增加
    colormap(jet);  % 也可用其他映射, 或自定义
    colorbar;
    xlabel('x_0');
    ylabel('y_0');
end
